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We examine the conditions for appearance of symmetry breaking bifurcation in damped and 
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phases exist between phases of symmetric normal and symmetric inverted oscillations. We find that 
symmetry broken solutions occupy a sufficiently smaller region of pendulum's parameter space in 
comparison to the statements made in earlier considerations [McDonald and Plischke, Phys. Rev. 
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I. INTRODUCTION 



Periodically forced and damped pendulum 



8 + 7$ + sin 8 = f cos uit 



(1) 



is one of the most important paradigms in the modern 
nonlinear science 1] . In solid state physics, the pendulum 
model well describes, for instance, the nonlinear dynam- 
ics of Josephson junctions of superconducting materials 

HI 

Fairly popular RSJ model of Josephson junction 0,13, 
which is equivalent to an underdamped pendulum, 7 <C 
1, demonstrates very rich nonlinear dynamics. Rota- 
tional states of the pendulum, (8) ^ (averaging over 
a period T = 2ir/w), correspond to a generation of dc 
voltage across the junction driven by ac current without 
dc component. Moreover, phase-locked states of the pen- 
dulum, (6) — n/m (n and m are integers), correspond to 
quantized values of dc voltage in the Josephson junction 
model. This effect, known as the inverse ac Josephson ef- 
fect 0, , has already found an application in the design 
of modern standard of one Volt y, H|- Chaotic vibra- 
tional and rotational motions of underdamped pendul um, 
7 <C 1 are also well-known in Josephson junctions |^llfj| 
(for review see |ll|). Surprisingly, the optimum operat- 
ing point for the voltage standards in (u), f )-parameter 
space is located near a region of chaos Therefore, 
knowledge of conditions for transition to chaos is impor- 
tant for the optimization of zero-bias voltage standard 
and related Josephson devices 0, Q • 

One of the most often observable roads to chaos is the 
period doubling scenario 0]. The pendulum 0} be- 
longs to a class of symmetric dynamical systems with 
the invariance under the transformations # — > — # and 
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t — > t + T/2. Therefore, symmetry-breaking bifurcation 
is a necessary precursor to period doubling 0, 0, [w| . 
This bifurcation describes a sharp transition from a sym- 
metric limit cycle satisfying 

8(t + T/2) = -0(t) + 2irk k is an integer (2) 

to symmetry broken limit cycles, for which this equality 
is invalid. As easily seen from QJ, a steady-state sym- 
metric solutions, 0(t), can have only odd harmonics of 
uj in their Fourier expansion series, together with a zero 
harmonic (8) that is equal to either or it. The solutions 
with (8) = are symmetric oscillations around the stable 
position 8 = 0, while (8) = ir corresponds to oscillations 
of inverted pendulum. In contrast, for symmetry broken 
solutions both even and odd harmonics are possible and 
(8) is some constant different from or w. Typically, 
the symmetry broken trajectories of underdamped pen- 
dulum occupy a very small range of parameters ojf just 
near transition to chaos [lol Il3l [po| . Detailed analyti- 
cal analysis of a transition from symmetric oscillations 
to symmetry broken oscillations near chaos border in the 
underdamped pendulum, 7 <c 1 and uj < 1, is presented 
in 

In spite of a widely spread belief that symmetry break- 
ing is always connected to a transition to chaos 1141 . this 
bifurcation is still possible for 7 -C 1 but uj > 1 |Tj|. In 
these conditions chaos is impossible and symmetry break- 
ing arises near transitions from normal to inverted states 
of the pendulum 01 • 

For a strong damping (7 > 1), neither chaotic [TsUIi^. 
nor rotational phase locked states ^(| can exist anymore. 
However, the presence of symmetry broken states of over- 
damped pendulum, that obviously are not related to a 
transition to chaos, has been briefly mentioned earlier in 
two pap ers devoted to dynamics of Josephson junctions 
p^L |20| . Moreover, it has been reported that these 
symmetry broken trajectories occupy a quite large part 
of the parameter space (oj, /). This interesting aspect of 
the pendulum dynamics at strong dissipation did not get 
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further attention so far, probably because of two main 
reasons. First, for the Josephson junctions pendulum's 
solutions at a strong damping are not physically interest- 
ing. Second, 9 and (9) do not correspond to any directly 
measurable physical variables in the junctions. 

However, it has been shown recently that another type 
of solid state microstructures - semiconductor superlat- 
tices subjected to a high-frequency electric field - can 
demonstrate very rich nonlinear dynamics similar to dy- 
namics of the Josephson junctions [H |H IHEa 
[27J ■ In particular, an analog of inverse ac Jose phs on effect 
has been predicted for the superlattices [23, l2£g • More- 
over, within some reasonable approximations nonlinear 
dynamics of an ac-driven semiconductor superlattice in 
the miniband transport regime is governed by a periodi- 
cally forced and damped pendulum |2(| |2?J • For the su- 
perlattices with realistic scattering constant [2^, [2t| , an 
effective damping in the corresponding pendulum model 
is not small: 7 > 1. Overdamped pendulum also arises 
in the models of lateral semiconductor superlattices [3(| . 

Importantly, in contrast to the case of Josephson junc- 
tions, a voltage across a superlattice is proportional to 
both the velocity, 9, and the coordinate, 9, of the pendu- 
lum (|TJ 26]. Therefore, even if rotations are impossible, 
(9) ^ 0, a dc voltage across the superlattice can still 
be generated due to contributions of symmetry-broken 
swinging oscillations with (9) ^ In (I = 0,1). For a 
strong damping, this is the only mechanism that can 
contribute to a rectification of THz signal in the semi- 
conductor superlattice |26|. 

The existence of physical situations, where the sym- 
metry breaking in pendulum at strong dissipation can 
be physically important, is the main motivation of our 
present work. Combining the analytical technique of 
truncated Fourier expansion [3l| with numerical simu- 
lations we find the conditions for symmetry breaking at 
strong damping. We describe a scenario of transition 
from symmetric to asymmetric oscillations at a strong 
damping. We found that symmetry broken (SB) oscil- 
lations form an intermediate stage between symmetric 
normal (N) and symmetric inverted (I) oscillations, i.e. 
N^SB^I^SB^N. Note that it is different from transi- 
tions N— >SB— >N with a large symmetry breaking phase 
even in overdamped case, which have been reported in 
|13| . We observed a relatively small symmetry breaking 
phase; with an increase of damping the ranges of driving 
amplitude and frequency, resulting in symmetry break- 
ing, decreases. Moreover, symmetry breaking does not 
exist in the overdamped case described by the first order 
differential equation. In this case only normal and in- 
verted oscillations survive. We presented simple formu- 
las providing a good approximation for bifurcation points 
between different types of symmetric (N— >I) and asym- 
metric (N— >SB, I— >SB) transitions in wide ranges of the 
frequency and the strength of alternating force. 

The organization of this papers is as follows. We start 
with a consideration of overdamped pendulum without 
inertial term (9 — 0). The analysis of this model appears 



to be the most simple and transparent. In the subsequent 
section irm we present analytic and numerical results on 
symmetry breaking bifurcations in the pendulum with in- 
ertial term. The final section of the paper is devoted to a 
summary and a brief discussion on applications in physics 
of semiconductors. The main text of our paper contains 
only most important results with sufficient explanations. 
Details of analysis and computing are presented in Ap- 
pendixes: Analytic analysis of bifurcations in first- and 
second order pendulum equations (Appendix^), numer- 
ical methods in use (Appendix 0) , proof of absence of 
hysteresis at symmetry breaking transition (AppendixQ 
and role of higher order harmonics at symmetry breaking 
(Appendix [D). 

II. FIRST ORDER PENDULUM EQUATION 

When damping is very strong we can neglect the sec- 
ond order derivative in Q and get the first order over- 
damped pendulum equation 

9 + (3 sin 9 = -q cos Lot (3) 

with (8=1/7 and r\ = f/"f. Numerical integration of 
this equation shows that only two types of stable periodic 
motion exists: Normal and inverted modes of oscillations. 
Solutions corresponding to these modes are all symmetric 
with (9) = and (9) = it for the normal and inverted 
mode solutions, respectively. Figure ^ shows a plot of 
regions of the two different modes in the (/, w)-parameter 
plane which exhibits a fan shaped pattern of alternating, 
disjoint areas of normal and inverted mode oscillations. 
A transition between the two modes is a sharp one and 
no intermittent states are in spite of rigorous attempts 
found. Especially, no symmetry broken regions exist. 

A simple analytic calculation verifies the absence of 
symmetry broken solutions and provides a good approx- 
imate condition for the transition. The analysis is done 
under the assumption that the solution is well described 
by the trial function 

9 = A + At cos(wi + at). (4) 

We begin by applying linear stability analysis to find the 
regions of stable symmetric motion. For details refer to 
Appendix JXJ. Straightforward analysis leads to a sta- 
bility condition 

(3 cos A J (At) > 0, (5) 

where Jo is the zeroth order Bessel function. If for the 
zeroth harmonic we have cosA < 0, then Jo(^4i) > 
must apply and vice versa. Both symmetric solutions 
lose their stability at the same At, which is a root of 
the Bessel function Jo. As the Bessel functions and the 
zeros of Jo will appear frequently, we adopt the shorthand 
notations 

Jo(^ (n) )-0, J { ^ = J k (A (n) ), (6) 
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FIG. 1: (Color) First order pendulum, with (3 — 1/2, regions of inverted and normal mode oscillations. White areas correspond 
to normal oscillations with (9) = and gray areas to inverted mode, {&) = 7T. No symmetry broken solutions are found. Solid 
(red) lines indicate the analytic prediction for mode boundaries, Eq. ©, which are plotted for n — 1 ... 6. 



where A^> is specifically the nth root of Jo and J£ n ' is 
Jfc evaluated at A^ . 

From the equation of motion, Eq. J2J, we can derive 
using the ansatz J2J relation 

J (A 1 )smA = Q. (7) 

From Eq. (J7J it follows that if A Q ^ 0, ir, then J (Ai) = 0. 
However, stability condition JSJ indicates that this case is 
not stable, as a perturbation does not decay. Therefore it 
is evident that asymmetric trajectories are not possible. 

Within the applicability of our first order approxima- 
tion, the only stable solutions are the ones with Aq = 
or 7r. Additionally, the pattern of normal and inverted 
mode solutions is revealed. Starting from a forcing r\ low 
enough, i.e. for which Jq(Ai) > 0, Eq. © implies that 
Aq = 0. When A\ crosses A^\ the sign must change 
in both factors, leading to Aq = iz. As 77 is further in- 
creased, A\ reaches A^> and the oscillations return to 
the normal mode. 

The boundaries of the normal and inverted oscillations 
can even be written in an explicit analytic form. From 
Eq. the expansion coefficient A\ in (@J can be found as 
a function of pendulum parameters. Substituting A\ — 
Ai n ) we get 

4l3 2 j[ n)2 +u 2 A^ 2 = v 2 . (8) 

Now we can compare our analytic results with numerical 
data. The boundaries of transitions between normal and 
inverted modes of oscillations, followed from Eq. 0, are 



superimposed on the numerical results in Fig. ^ We see 
that the the equality Jq(Ai) = and the prediction of 
Eq. JSJl hold well at the transition lines with reasonable 
accuracy everywhere, except in the case of low frequency. 
In the limiting case of very low frequency drive, the trial 
solution in the simplest form Q can become invalid and 
effects of higher harmonics should be taken into account. 

III. SYMMETRY BREAKING IN THE 
PENDULUM AT STRONG DISSIPATION 

A. Numerical integration 

Now we return to the pendulum with inertia, Eq. Q • 
We start with a review of our numerical results; for a de- 
tailed description of numerical methods see Appendix [H] 
The study of (to, /)~parameter plane reveals a fan shaped 
structure of regions of normal and inverted mode oscil- 
lations, similar to that of the first order pendulum. In 
addition to the alternating pattern of the symmetric so- 
lutions, narrow regions of symmetry broken solution have 
now emerged on the interfaces between the regions. For 
the solutions in these regions we have observed nonzero 
values of zero harmonic components along with nonzero 
even harmonics, which are requisites of a symmetry- 
broken solution. The different modes on the parameter 
plane are depicted in Fig. [21 where we have chosen 7 = 2 
for the damping coefficient. Symmetry broken regions 
appear around both transitions, from normal to inverted 
and back, thus separating the normal and inverted mode 



4 




FIG. 2: (Color) Second order pendulum equation, regions of inverted and normal state oscillations. Sandwiched between them 
are the regions of symmetry breaking. Solid (red) lines indicate analytic prediction of Eqs. (11 2t and ill IH . 



regions completely. Regions are not wide, in contrast to 
what was reported in |13( |4J| . Increasing the damping 
will cause them to get more narrow. A closer examina- 
tion is shown in Fig. where we have fixed ui = 0.2 
and 7 = 3. Harmonic components Aq and A\ are plotted 
against the drive amplitude /, and the sequence of modes 
of solutions is made more clear. Normal mode bifurcates 
into asymmetric solutions when A\ is close to a root of 
the Bessel function Jq(Ai). This symmetry broken state 
persists only for a short interval of variable /, ending in 
inverted mode, when Aq has reached n. Small deviations 
from approximately linear response of A\ to a change of 
/ can be seen at the transitions from normal to inverted 
mode. 

Close-ups of the symmetry-broken range provide more 
insight. The following two examples are depicted in 
Fig. The first case we consider is of low frequency 
drive, u> = 0.2 and 7 = 3. Symmetry breaking starts for 
moderate forcing, / = 1.8413 and with Ai slightly larger 
than A^\ the first root of Jo. The curve of Aq of sta- 
ble symmetric solution bifurcates from zero into a curve 
corresponding to a stable asymmetric solution. This is 
accompanied with the emergence of even harmonics in 
the Fourier spectrum and a sudden stop in the increase 
of A\ . Further increases of / cause A\ to decrease, while 
Ai traces a sugar-loaf shaped curve. It should be noted, 
that the curve of the second harmonic reaches quite high 
values, Ai ~ 0.6 at the maximum. The zero harmonic 
varies continuously and monotonously from zero to ir. 
We should note that actually two different stable solu- 
tions exist in the symmetry breaking range, that are re- 
lated to each other by the symmetry. Initial values of the 



system determine which branch the solution converges to 
after the initial transient. Symmetry breaking ends in an 
inverted state at / = 1.8773 with A\ just less than A^\ 
In the second case we examine, we set uj = 4 and 7 = 3. 
To achieve symmetry breaking, one needs high values of 
/. In this case we have / = 47.46. Again, A\ has just 
crossed A^> when the symmetric mode loses stability. 
Now the decrease of A\ in the symmetry broken region 
is less dramatic than in the low frequency case. Further- 
more, the second harmonic component is almost negligi- 
ble. In other respects, this case doesn't differ from the 
previous example. 

B. Analytic analysis 

In this subsection we present explicit analytic form of 
conditions for the bifurcation to occur and provide a pic- 
ture of the scenario of symmetry breaking that is con- 
sistent with the numerical results. Here, as in Sec. [HJ 
we assume that oscillations of the pendulum can be de- 
scribed by the trial function in the form For details 
and discussion on the range of validity of our analytic 
results, see Appendix lAl 

Following the method of [Tol l3l| , we linearize Eq. 
and, with appropriate approximations, put it into the 
form of the Mathieu equation. Using the existing knowl- 
edge of the solutions of Mathieu equation, we derive the 
following conditions for bifurcations from symmetric os- 
cillations 

J (A 1 ) + (-l) fc 2 2 J 2 (A!) 2 =0, (9) 
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FIG. 3: (Color) The amplitudes of the two first harmonic components, Ao (blue) and A\ (red) plotted against /. Ai shows 
nearly linear dependence on parameter /. Small deviations from this are seen near the transitions between normal and inverted 
states, along with a Ao component that is not equal to 0,tt. To guide the eye, dashed lines connecting the points Ai(f) = A^ 
on the vertical axis to the corresponding / on the horizontal axis are drawn. The fixed pendulum parameters are 7 = 3 and 
lo = 0.2. 



where k = corresponds to a normal mode, while k = 1 
marks an inverted mode. Notice, that Eq. © is the same 
condition for loss of stability as the one for the first order 
pendulum, Jq(Ai) = 0, but now with an additional term. 
This term only becomes significant when Ai is close to 
a root of the Bessel function J . As a consequence, now 
there are two separate bifurcation points corresponding 
to the two different symmetric modes. These are located 
rather close to the roots of Jo. Clearly, the dynamics of 
the second and the first order pendula are similar away 
from the regions of transition, justifying the correspon- 
dence to the analysis of the previous section. 

The explicit solution of Eq. (0 in terms of amplitude 
A\ can now be found. Owing to the fact that the latter 
term in Eq. @ is small, we get for the critical values of 



A N (n) 

X 



= A^± 



2. J. 



j{">(f + 4 w 2) Tx («) 



(n) 



(10) 
(11) 



Here A^ n > and Ai^ n ' are respectively the critical ampli- 
tudes of normal and inverted modes for which symmetric 
oscillations become unstable. A small asymmetry with 
respect to A^ in the critical A\ due to the term \^ can 
be neglected within the applicability of Eq. H10[l without 
loss of qualitative agreement. This is implicitly assumed 
hereafter in this section. 

In transitions from normal to the next inverted mode, 
critical value of amplitude A\ for the inverted mode is 



less than that of the normal mode: A\- n) < A N ( n \ The 
reverse applies for transitions from an inverted phase to 
the next normal phase region, i.e. A N (n) < Aj( n \ There- 
fore Ai must decrease as the boundary of the first sym- 
metry broken region is crossed, and the nth region starts 
when Ai > A^ and ends when A\ < A {n \ Such kind 
of behavior of A\ is in a good agreement with numerical 
results shown in Fig.0] Moreover, following Eq. (|10fl the 
difference between An'"' and A^ n ' decreases with an in- 
crease of the frequency to. This explains why the range of 
A\ within the symmetry breaking region is wider in the 
low frequency case in comparison with the case of high 
frequency drive (cf subplots (a) and (b) in Fig.^J. 

We turn now to finding of surfaces in the parameter 
space (7, bj, /), on which the symmetry breaking bifurca- 
tion happens. First, we substitute the trial solution (JIJ 
in the equation of motion Q and in the first harmonic 
approximation derive the dependence of amplitude A\ 
on the pendulum parameters 7, ui and /. Second, sub- 
stituting critical values of A\ from Eq. (|1(J|) and taking 
Aq = kir, we have 



f 2 = (jujA,) 2 + (lu 2 A 1 ±2J 1 (A 1 )Y 



At 
Ax 



In) 
» 



(12) 

for bifurcations from inverted and normal modes, respec- 
tively. Taking the difference of these two critical / does 
indicate that one symmetric mode loses stability before 
the other gains it. This eliminates a possibility of hystere- 
sis around the symmetry breaking transition and shows 
that both symmetric phases cannot be stable for same 
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(a) Low frequency and weak drive, 7 = 3 and 
u) = 0.2. 



(b) High frequency and strong drive, 7 = 3 and 
a; = 4. 



FIG. 4: (Color) The amplitudes of the two first harmonic components as functions of the drive amplitude /. Notice the 
relatively large amplitude of the second harmonic component in the low frequency case. 



pendulum parameters. For more detailed discussion on 
this subject, see Appendix IO 

Our main analytic results, Eqs. i|12[) and (|10l) . are ex- 
pected to be valid for uo > 1. Even for frequencies less 
than unity, these results may be applied, provided that / 
is sufficiently large. Increasing the damping does some- 
what improve results, and for low frequency overdamped 
case, good qualitative picture can be still provided even 
for low frequencies, to ~ 0.1 (see Appendix |A"|) . The 
effects of third and higher order harmonics in the trial 
solution (@J are discussed in Appendix iDl 

On the other hand, Eq. I|12|) shows remarkable agree- 
ment with the numerical results for a very wide range of 
parameters. In Fig. [5] analytic prediction for the bound- 
aries of symmetry breaking regions are shown as solid 
(red) lines. Additionally, in Table [I] we have tabulated 
points of symmetry breaking bifurcations in / variable 
for a selected set of quite different values of ui and 7. 
Boundaries of regions of stable symmetric solutions can 
be found even in low damping and moderately low fre- 
quency case, provided that / is sufficiently large. This 



suggests that as damping or frequency are lowered from 
sufficiently high values, the symmetry broken regions 
evolve into the regions of complex and chaotic dynam- 
ics with the alternating normal and inverted phase mo- 
tion separating them. Evidently, our analytic results pro- 
vide a reasonable approximation even in parameter space 
regimes beyond only overdamped case. 



Now we can summarize the main results of present 
section. Taking the alternating force strength, /, as the 
control parameter, we can report the following behavior: 
A stable inverted or normal mode becomes unstable at 
an amplitude of stationary oscillations, A%, that is just 
greater than A^ and symmetry breaking starts. Increas- 
ing / further causes Ai to decrease, as the bifurcation to 
the next mode should occur when A\ is less than A^ n '. 
The critical values of the A\ are given by Eq. I|1(J|I with 
the critical parameters by Eq. (|12|l . These critical points 
approach a common value, A^ n \ as damping is increased. 
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Parameters Numerical Analytic Errors (%) 
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n 


h 


h 


h 


h 


eMax 


ewidth 


2 


1 


1 


5.19964 


5.76731 


5.22416 


5.7492 


0.48 


7.6 


3 


2 


2 


39.4625 


40.1666 


39.4589 


40.1622 


0.011 


0.11 


0.05 


2 


1 


8.76666 


10.51386 


8.79088 


10.5046 


0.28 


2.0 


0.5 


1 


2 


5.61680 


6.79113 


5.60002 


6.76786 


0.35 


0.56 


0.5 


0.5 


3 


2.73651 


3.53534 


2.71429 


3.45881 


2.2 


6.8 



"Region between the critical / involves in addition to symmetry 
broken solutions also complex dynamics including chaos. 

TABLE I: Locations and widths of selected symmetry breaking intervals. Here {^,uj,n) are the parameters of the system with 
n being the number of the transition. The start and end points of the symmetry broken range are, respectively, /i and /2. 
In the errors' column, eMax is the larger of the relative errors in the start or end points of the symmetry breaking interval. 
Additionally the error in the width of the symmetry broken interval ewidth is given. Analytically given values for the start and 
end points are obtained using Eq. 112H along with Eq. 1101 . Notice the good agreement between the analytic and numerical 
result even for very low damping and moderate frequencies. 



IV. CONCLUSION 

In summary, we confirmed the existence of symmetry 
breaking phenomenon in the strongly damped pendulum. 
The symmetry breaking states form a necessary stage in 
transitions from the normal to the inverted and from the 
inverted to the normal states of the overdamped pendu- 
lum. The scenario of transitions from the symmetric to 
the symmetry broken states, and vice versa, appears to 
be similar to the case of underdamped pendulum without 
chaos O ur results also demonstrate that an inertial 
term in pendulum is a requisite of symmetry breaking. 
On a more hypothetical note, it could then be possible to 
achieve symmetry breaking similar to what it described 
here, if higher even derivatives or/and additional nonlin- 
ear damping terms are included. 

These results can be directly applied to the studies of 
novel mechanism of THz radiation rectification in semi- 
conductor superlattices [2f| because this system has un- 
derlying pendulum dynamics pft l3(il | . Note that in some 
cases pendulum-like equation of the third order j32, 
and a pendulum equation with a nonlinear damping term 
supporting symmetry of the problem .34j naturally arise 
in different models of semiconductor superlattices. We 
speculate that swinging symmetry broken oscillations in 
these superlattice models also correspond to a sponta- 
neous generation of dc voltage or dc current. Related 
detailed consideration will be published elsewhere. It 
is worth to notice also the importance of the symme- 
try breaking bifurcations in the physics of ac-driven bulk 
semiconductors with different types of nonlinearity |35| . 
Our results can be also useful in the development of this 
area of semiconductor nonlinear dynamics. 
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APPENDIX A: DETAILS OF ANALYTIC 
ANALYSIS 

1. First order pendulum 

We start with a derivation of stability condition J5j|. 
Consider a small perturbation 8(t) to a solution of Eq.©. 
The evolution of the perturbation is governed by the lin- 
earized equation 

5 + f3cos9 5 = 0. (Al) 
Its formal solution is 

S = exp (-ft J dt cos 6^ . (A2) 
Using Eq. (0} and well known identities 

oo 

sin(acos</>) = 2 ^(-1)" J 2 „+i(a) cos ((2n + l)</>) , 

(A3) 

OO 

cos(acos0) = J (a) +2 VV-l) n J 2n (a) cos(2n0), (A4) 

n=l 

the solution S(t) is put into the form 

S = exp (—f3 cos A J (Ai)t + periodic terms) . (A5) 

A perturbation to a stable solution must be damped, 
which leads to the stability condition JSJl. 

Now we turn to the derivation of Eq. JSJ • We substi- 
tute the ansatz into Eq. ©. Identities (|A3|> and IIA4I) . 
truncated to include terms up to the first harmonic, are 
now employed to expand the sine term. Finally, in the 
resulting expression equating the trigonometric functions 
separately along with the zero harmonic, the following set 
of equations is obtained. 

— 2/3cos A Ji(Ai) = r]cosai (A6a) 
loAi — ?7sinQ!i (A6b) 
Jo(Ai)sinA = (A6c) 



8 



Eliminating the phase between Eqs. (|A6a|) and l|A6b|l we 

get 

4/3 2 cos 2 yloJi(Ai) 2 +a; 2 ^ = v 2 . (A7) 
Eq. JHJl follows by substituting A\ = and A = 0, ir. 

2. Second order pendulum 

a. Derivation of Eqs. 0) and 11(A) 

Linearization of Eq. Q yields the equation for the evo- 
lution of an infinitesimal perturbation 



<5 + 7(5 + 008 6* S = 0. 



(A8) 



Change of variables S = exp(—jt/2) £ removes the first 
order derivative term. The trial function (@J is again 
used to substitute 6 and the cosine is expanded up to 
second order, as the first order term vanishes for sym- 
metric mode solutions. Finally, after rescaling the time 
t = ujt + ai, Eq. I|A8J) is put into the form 



(-1)* 



Jo(^) - {-l) k \ - 2J 2 {A 1 ) cos2r 



e=o. 

Overdot stands for differentiation with respect to r and 
we use Ao = kn. It is easy to see that Eq. (IA9(I has the 
form of the Mathieu equation, 



y + (a — 2q cos 2r)y = 0, 
when we make the identifications 

_ (-l) k 4Jo(A 1 )- 1 2 
4a, 2 

q = { -^-J2(A 1 ). 



(All) 
(A12) 



Following Floquet theorem, a solution of Mathieu equa- 
tion can be presented in the form y = P(r)exp(wr), 
where P is a periodic function and v is the Mathieu 
characteristic exponent. Stability of the symmetric os- 
cillations is now determined by the condition 



J_ 

' 2uj 



SJm v < 0, 



(A13) 



which holds for stable solutions. Therefore, a symmetric 
solution loses stability at v — —17/20,. Real part of v is 
zero. 

We now turn to finding the points of bifurcation. We 
can make use of the expression for a as a function of q 
and v 361 as 



v 2 + 



2{v 2 - 1) 



(A14) 



which we truncate to include only the two first terms. 
Substituting the expressions for ^, a and q into (|A14|) . we 



obtain Eq. which is the condition for a bifurcation 
from a symmetric solution. Note that the truncation of 
series (|A14I) is valid provided that q is small or v 2 is large. 
These conditions are usually satisfied, independent of lu, 
because higher order terms are all proportional to 1/lo 2 . 
Large damping is then needed for small lo in order to 
have the series converge fast enough for the truncation 
to apply. 

Equation © can be approximately solved. We substi- 
tute A\ — A^ + £ and make a series expansion up to 
the first order in e. 



2(-l)*J. 



k Any 



( 7 2 + 4c 2 )J 1 ( " ) + (-l)V" ) 



= 0, 



(A 15) 

where X (n) is defined in Eq. (JTT1)- The validity of this 
approach is due to the fact that the second term in Eq. @ 
is small for wide range of parameters. It is obvious that 
7 2 + lo 2 3> 1 is sufficient. However, this can be even 
slightly relaxed. Symmetry breaking is expected when 
A\ is near a root of Jo- Now, due to known properties 
of Bessel functions, we may infer that if A\ is sufficiently 
large and near a root of Jo, it is also close to a root of 
J%. Thus, the discussed second term may be small even 
for low values of 7 and lo. 

Resulting equation i|A15|) is trivial to solve for e. So- 
lution yields 



e = 



2(-l 



i k J, 



J 1 (n) (7 2 + 4o, 2 )-(-l) 



fcy(n) 



(A16) 



(A10) from which Eq. llOf) immediately follows. Note that the 



X 



(n) 



in the denominator results in a slight asymmetry 
of the bifurcation points with respect to A^ n \ Within 
the applicability of our results, 

x («) i s 

small enough to 

be neglected, when a more tractable form of Eq. (|A16|) 
is needed. Some of the quantitative agreement is lost 
in this approximation, but qualitative the result is the 
same. This is due to the fact that in the derivation of 
Eq. (1 A 15(1 we already assumed, that j 2 +oj 2 is sufficiently 
large, or that J2{Ai) is small. 



6. Derivation of Eq. 112ft 

Substituting the trial solution 8 = A + A\ cos(ojt + 
«i) into the Eq. JIJ, expanding the sin6> term using the 
second order truncations of identities i|A4|l and i|A3|l and 
finally equating the harmonic terms, the following set of 
equations is obtained 

sin^o Jo(Ai) = 0, (A17a) 
2cosA Ji(Ax) - lj 2 Ai = /cosai, (A17b) 
— juiAi = /sinai. (Al7c) 



Eliminating phase between Eqs. i|A17b(l and l|A17c|) we 

get 

f 2 = (7WA1) 2 + (w 2 ^i - 2 cos ^0 Ji{Ai)) 2 . (A18) 
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Replacing Aq and Ai with values corresponding to the 
bifurcations of normal and inverted modes, we end up 
with an explicit expression for critical points, Eq. 1)12(1 . 



APPENDIX B: ON NUMERICAL METHODS 

In all numerical calculations, standard double precision 
floating point arithmetic has been used. For the numeri- 
cal integration of system , the equation is put into the 
form 



71 = 02, 

9 2 = — 7#2 — sin#i + / cosivt, 



(Bl) 



where 6\ = 9. Initial values 9 k {0) = 0, k = 1, 2, were used 
in all computations. The 4/5th order embedded Runge- 
Kutta routine with Cash-Karp coefficients [37j is used in 
majority of computations in this paper. For the compar- 
ison of analytic and numerical results on the symmetry 
breaking regions, Table |U 7/8th order Runge-Kutta with 
Prince-Dormand coefficients |3Sj was employed. In vari- 
able step-size control, a relative accuracy of around 10~ 10 
was typically used. 

Harmonic components were extracted from the numer- 
ical solution by introducing additional differential equa- 
tions of the form 

£o = £fc = cos ku>t, (k=9smku>t, (B2) 

for k > 0. Same error bounds and step functions were 
used as for the variables 9\ and 82- Coefficients in the 
Fourier series expansion 



6 = a + (afc cos ku>t + b k sin kut) 



k=l 



have been evaluated using 
1 



ao = jj(fo(to+T)-eo(to)), 
a k = |(&(to+T)-&(*o)), 
b k = |(C*(*o+T)-C*(to)). 



(B3) 

(B4) 
(B5) 
(B6) 



This simply amounts to calculating the projection of the 
solution onto orthogonal base of trigonometric functions. 
Finally, cosine harmonic coefficients of the series 



9 = Aq + A k COs(kujt + Qffc), 



k=l 



are computed using 
Ak = \l 

-Oik = 



4 + bl 




arctan 1 






>:: 


arctan I 








bk t 




M 2' 






(B7) 



(B8) 



(B9) 



Initial transient phase was eliminated by first performing 
a period search on the solution. A simple integration in 
steps of one drive cycle, T — 2it/uj, was made and after 
each step the current values of 9 and 9 were recorded and 
compared against the values of previous steps. A period 
of nT is assumed, if 

\9 k (t) - 9 k {t - nT)\ < e{\e k (t)\ + \e k (t-nT)\), 



applies for k — 1,2 (or just k — 1 for the first order 
pendulum). Tolerance is determined by the parameter 
e. Typically, a value of e = 10~ 9 was used. When near 
bifurcation points, lower values were needed, e = 10~ n 
with correspondingly higher integrator accuracy require- 
ments. 

To determine numerically the bifurcation points re- 
ported in Table [I] a finite time maximal Lyapunov ex- 
ponent was used |3fll l4fT | . Because one requires possibly 
a very long integration time to reach the stable limit cycle 
in the vicinity of bifurcation points, it was not deemed 
feasible to find the actual values of the harmonic com- 
ponents. Rather, Lyapunov exponent as a measure of 
the rate of convergence of a small perturbation is good 
indicator of bifurcation points, as it tends to zero as the 
critical point is approached. 



APPENDIX C: ABSENCE OF HYSTERESIS 



In Section lill Bl it was stated that there is no hystere- 
sis around symmetry breaking transition, because one 
symmetric mode loses stability before another gains it. 
However, actually situation is more complicated because 
regions of stability of symmetric normal and inverted 
modes do overlap, if one considers just A\ as the vari- 
able. In this Appendix we will show that there is no 
contradiction between overlapping of regions of stability 
of the symmetric modes in (j4o> ^4i)-plane and absence of 
such overlapping in (/, w)-plane. 

To begin with we want to demonstrate that the stabil- 
ity regions in units A\ really can overlap. If the asym- 
metry of the bifurcation points with respect to A^ in 
Eq. (|A16II is neglected, the critical values of A\ for sta- 
bility are given by 



Al (n) 



= A^ ± e A 



T (nV 



7 2 + Au> 2 jv 



(n) 



(CI) 



(C2) 



Clearly, ea < for even and ea > for odd n. As odd n 
correspond to transitions from normal to inverted mode 
(N— >I) and the other way around (I— >N) for even n, the 
stability regions overlap. 
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Now we consider the difference between criticai points 
of stability in variable /. We substitute (|(J1|) in (1 A 18(1 
and expand in ea up to first order, then substitute k 
and Ai with values corresponding to inverted and normal 
modes. Taking the difference of these values gives 

A/ 2 = //"^-/n'") 2 (C3) 

+8J i 1 n) cj 2 A^ + SJ^J^sa. (C4) 

We next drop the last term, as it is small under the as- 
sumptions we have made earlier in derivation of Eq. I|(J1[1 . 
Overlapping can occur when A/ 2 is negative for n = 

1,3,... and when A/ 2 is positive for n = 2,4, For n 

being odd this condition results in 

7 2 + 4w 2 2 1 ' 1 j 

which holds if regions of stability overlap. Estimating left 



hand side upwards, (j 2 + Lo 2 )j(j 2 + 4uj 2 ) < 1, we have 

> J\ , which is false. Analogously the case of 
even values of n can be considered. These considerations 
confirm that both symmetric modes cannot be stable si- 
multaneously for same values of parameter /. 



APPENDIX D: APPLICABILITY OF THE FIRST 
HARMONIC TRIAL FUNCTION 

The analytic results derived herein all rely upon the 
assumption that 9 = A + A 1 cos(u>t + a\) is a good 
approximation of a correct T-periodic solution or at 
least captures the essential features of it. Numeri- 
cal integration of systems and © indicate that at 
least when to > 1 applies the higher harmonic com- 
ponents are insignificant. For analytic analysis we in- 
cluded the third harmonic term into the trial function, 
8 = krr + Ai cos(ujt + ai) + A3 cos(3ojt + as) . Substituting 
it in motion equation Q we get 



6(-l) fc ^J (A 1 ) ^ (A 3 ) + (-l) l 2u, ^-^MA^MArf + ( 7 2 + 9uj 2 ) J (A 3 ) 2 J 3 (A^ 

3w 2 ( 7 2 + 9lj 2 ) 



(Dl) 



r 



where I is an integer determined by conditions not es- 
sential to current analysis. Eq. (|D1|I vanishes as 1/lu 2 
for high frequencies, in support of the numerical result. 
It is noteworthy that A±, and consequently /, are cou- 
pled to A3 only via the Bessel functions in the numer- 
ator. Thus, high values of / and A\ cannot increase 
significantly the third harmonic. Strong external drive 
therefore makes the relative contribution of the higher 
harmonics smaller. Clearly, problems can arise when the 
frequency is very low. Numerical integration also shows, 
that higher harmonics become more significant and even 
Eq. (|D1(I breaks down. Although individual harmonic 
components may have relatively low values, we have ob- 
served that the spectrum as a whole can be quite spread. 
Therefore, in this limiting case an inclusion of the first 
two or three harmonics into the trial function may not 
be enough. 

Now we discuss the influence of damping. Surprisingly, 
we have found that in general an increase of damping 
does not always suppress the higher harmonics. From 
Eq. HD1I) we have 

a 3 « — vjo(^ 3 ) 2 j 3 (^i) 2 - mWj^m 2 (D2) 

W7 



for very large damping. A3 can be large even for large 7, 
if frequency is low. The Bessel functions become signifi- 
cant and strongly nonlinear response to an increase of 7 
is expected. Numerical integration confirms this. Rela- 
tive amplitude of the third harmonic with respect to A\ 
fluctuates with an increasing amplitude as 7 is increased. 
Eventually these fluctuation will die out. There is an up- 
per limit for 7, approximately given by (|12fl when n = 1, 
for which a symmetry breaking bifurcation can happen. 
We have observed numerically that after this final critical 
7 is crossed, fluctuation stop and the higher harmonics 
slowly fade away. 

In summary, the first harmonic approximation of the 
solution is definitely valid when ui j> 1. Increasing damp- 
ing alone will not make the truncation more applicable in 
the regime, where interesting phenomena are expected. 
It will improve the accuracy of some results presented 
in this paper, which rely on assumption that lo 2 + 7 2 
is large. Specifically, overdamped case is well described 
even for lo ~ 0.1, albeit without impressive quantitative 
agreement. High values of parameter / result in higher 
relative proportion of the first harmonic component in 
the solution, thus improving the approximation. 
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